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Abstract 

This article presents a new distance for measuring shape dissimilarity 
between objects. Recent publications introduced the use of eigenvalues of 
the Laplace operator as compact shape descriptors. Here, we revisit the 
eigenvalues to define a proper distance, called Weighted Spectral Distance 
(WESD), for quantifying shape dissimilarity. The definition of WESD is 
derived through analysing the heat-trace. This analysis provides the pro- 
posed distance an intuitive meaning and mathematically links it to the 
intrinsic geometry of objects. We analyse the resulting distance defini- 
tion, present and prove its important theoretical properties. Some of 
these properties include: i) WESD is defined over the entire sequence of 
eigenvalues yet it is guaranteed to converge, ii) it is a pseudometric, iii) 
it is accurately approximated with a finite number of eigenvalues, and iv) 
it can be mapped to the [0, 1) interval. Lastly, experiments conducted on 
synthetic and real objects are presented. These experiments highlight the 
practical benefits of WESD for applications in vision and medical image 
analysis. 

1 Introduction 

Quantifying shape differences between objects is an important task for various 
areas in computer science, medical imaging and engineering. In manufactur- 
ing, for example, one may wish to characterize the difference in shape of two 
fabricated tools. In radiology, a doctor frequently diagnoses a disease based 
on anatomical and pathological shape changes over time. In computer vision, 
discriminative shape models are used for automated object recognition, [T7|[44]. 

In order to define measurements of shape dissimilarity, scientists rely on 
descriptors of objects that capture information on their geometry 44 . These 
descriptors can be in the form of parametrized models (e.g. point clouds, sur- 
face patches, space curves, medial axis transforms) or in the form of geometric 
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properties (e.g. volume, surface area to volume ratio, curvature maps). Once 
a descriptor is formulated the distance between two shapes can be defined as 
the difference between the associated descriptors. The exact definition of the 
distance however, is a critical issue. In order to define an intuitive and theoret- 
ically sound distance, one should ensure that it takes into account the nature 
of the descriptor. For instance, the descriptor might be an infinite sequence of 
positive values, in which case we should be careful not to define a distance that 
diverges for every non-identical pair of shapes. 

Shape descriptors based on the eigensystems of Laplace and Laplace-Beltrami 
operators, called spectral signatures, have recently gained popularity in compu- 
tational shape analysis [8j[23j[33j[34j[35j[38] . These descriptors leverage the fact 
that the eigenvalues and the eigenfunctions of Laplace operators contain infor- 
mation on the intrinsic geometry of objects [iT 20,4l]. A visual analogy useful 
for an intuitive understanding is to think of an object (e.g. in 2D) as the mem- 
brane of a drum. In this case the eigenvalues correspond to the fundamental 
frequencies of vibration of the membrane during percussion, and the eigenfunc- 
tions correspond to its fundamental patterns of vibration. Both the eigenvalues 
and the eigenfunctions depend on the shape of the drum head and thus can be 
used as shape descriptors for the object. 

Despite recent progress by [8j|23j|33[|34j[35j[38] , designing meaningful shape 
distances based on spectral signatures remains challenging. Difficulties arise 
from the nature of the eigensystems. The eigenfunctions of a shape mostly pro- 
vide localized information on the geometry of small neighborhoods. Aggregating 
such local information into an overall shape dissimilarity measure is non-trivial. 
On the other hand, the eigenvalues provide information about the overall shape, 
so they are ideal for defining global distances. However, they form a diverging 
sequence making it difficult to define a theoretically sound metric. Here, we 
tackle this latter problem and propose a new shape distance based on the eigen- 
values, which is technically sound, intuitive and practically useful. 

In the remainder of this section, we first review in further detail the liter- 
ature on spectral signatures and shape distances related to eigenfunctions and 
eigenvalues. Then, we provide a brief overview of our new shape distance. 



1 . 1 Eigenfunctions 

The eigenfunctions of an object constitute an infinite set of functions. Each 
function depends on the shape of the object and is different than the rest of the 
set. Figure [l] illustrates this for two example objects where a few eigenfunctions 
are shown. The values these functions attain at each point capture the local ge- 
ometry around the point, i.e. of its neighborhood. Inspired from this geometric 
information, methods define local shape signatures [18^,23^^35^,38^ for each point 
on an object by evaluating a subset of eigenfunctions at that specific location. 
Global shape distances are then defined using such local signatures. Such dis- 
tance definitions rely on correspondences. These correspondences should hold 
both in terms of points and the subset of eigenfunctions used in the signatures, 
a condition hard to satisfy in practice ^18^. Explicitly searching for such corre- 
spondences leads to expensive algorithms [5j[9j[T8j[28^^30^ . On the other hand, 
computing distances between distributions of local signatures obtained by ag- 
gregating all the points, as in [22{|35[ |38^ , might implicitly construct false corre- 
spondences. In summary, defining a global distance based on local signatures is 
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Figure 1: Starfish and tarantula. The objects represented as binary maps are 
shown on the left, followed by the 1st, 2nd, 5th, 20th, and 100th eigenfunction. 
The values increase from blue (negative) to red (positive) with green being zero. 



not an easy task. 

Instead of extracting local information from an eigenfunction, one can also 
think of capturing its global pattern by looking at regions where its values are 
all positive or all negative. Such regions are called nodal domains. Different 
eigenfunctions induce different patterns and, in turn, have different number 
of nodal domains, called nodal counts 



11 . For a given object, the ordered 



sequence of nodal counts contain information on its overall geometry 13 



Inspired by these observations, authors in 21 used this sequence as a global 
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shape signature. They further defined the associated shape distance between 
two objects as the Euclidean norm of the vector difference between their nodal 
count sequences. However, it is not intuitively clear what the nodal counts 
represent. Furthermore, the entire sequence is diverging so that, in practice, 
one first chooses a finite subset and then computes the distance for that subset. 
These difficulties make it hard to define an intuitive and sound shape distance 
based on nodal counts. 



1.2 Eigenvalues 

Signatures based on eigenvalues, on the other hand, have a clearer geometric 
interpretation. The set of eigenvalues contains information on the overall ge- 
ometry of the object. Specifically, the ordered sequence is analytically related 



to the intrinsic geometry by the heat-trace^ [27 31 32,37,40. Hence, more 
intuitive distances can be constructed using the eigenvalues. However, simi- 
lar to the sequence of nodal counts, the eigenvalue sequence is also divergent. 
This makes the distance definition theoretically challenging. Inspired by the 



sequence's link to the geometry. Renter et al. in 33 , used the smallest N 
eigenvalues as a shape signature, called shape-DNA. As the associated shape 
distance, the authors proposed the Euclidean norm of the vector difference be- 
tween the shape-DNAs of objects. Although this is a very good first attempt 
the divergent nature of eigenvalue sequence results in important theoretical lim- 
itations for this distance, as also pointed out in |28 . The main problems are 
i) defining a distance on the entire sequence does not yield a proper metric, ii) 
the differences between the higher components of two sequences dominate the 
final distance value, even though these components do not necessarily provide 
more information on the geometry, and iii) the distance value is sensitive to the 
choice of the signature size N. These theoretical problems also cause practical 
drawbacks as we demonstrate later. 
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This article proposes a new shape distance, cahed Weighted Spectral Dis- 
tance (WESD), using the sequence of eigenvalues of the Laplace operator. We 
derive WESD from the functional relationship between the eigenvalues and the 
geometric invariants as given by the heat-trace. This derivation provides WESD 
a clear geometric intuition as a shape distance. It also links WESD to the dis- 
tance defined by Renter et al. in fSSi as well as to the local signature defined 
in 35 . The resulting formulation of WESD differs from other previously pro- 
posed scores based on eigenvalues, whether in shape analysis or other fields |l9], 
both in its formulation and in the fact that it is defined over the entire sequence. 
This latter point, as we will show later, alleviates the critical importance of the 
choice of the signature. We furthermore analyse and prove theoretical properties 
of WESD showing that it does not share some of the fundamental problems the 
distance proposed in |j33j has. Specifically, we prove that WESD: i) converges 
despite the fact that it is defined over the entire eigenvalue sequence, ii) can be 
mapped to the [0, 1) interval, iii) is accurately approximated with a finite num- 
ber of eigenvalues and the truncation error has an analytical upper bound and 
iv) is a pseudometric. These theoretical properties also yield important prac- 
tical advantages such as being less sensitive to the signature size (truncation 
parameter) providing a principled way of choosing this parameter, providing 
more stable low-dimensional shape embedding and simplicity in combining with 
other distances as WESD can be normalised. Applying to synthetic and real 
objects, we further demonstrate the benefits of WESD in comparison to the 
other eigenvalue-based distance defined in |33j . 

The remainder of this article is structured as follows. Section |2] presents 
a brief overview of the Laplace operator, the eigenvalue sequence and its role 
in shape analysis. In Section [3] we define WESD and derive its theoretical 
properties. Section |4] presents an extensive set of experimental analysis on 2D 
objects extracted from synthetic binary maps, shape-based retrieval results for 



3D objects using the SHREC dataset 25 , low dimensional embeddings of real 



3D data such as subcortical structures in brain scans and 4D analysis of binary 
maps extracted from cardiac images. 



2 Spectrum of Laplace Operator 

This section provides a brief background on the Laplace operator, its eigenvalue 
sequence, called spectrum^ and its role in shape analysis. We first relate an ob- 
ject's intrinsic geometry to the spectrum of the corresponding Laplace operator. 



and 33 



We then provide some details on the previously proposed shape-DNA 33 and 
discuss the associated issues. For further details we refer the reader to ll|32pQ 



We denote an object as a closed bounded domain ^7 C M with piecewise 
smooth boundaries. In the case of binary maps, Q would correspond to the 
foreground representing the object. For a given the Laplace operator on this 
object is defined with respect to a twice different iable real- valued function / as 



^^f = l^d^f^ ^^^^ 

i=l * 

where x = [xi, ...^Xd] is the spatial coordinate. The eigenvalues and the eigen- 
functions of are defined as the solutions of the Helmholtz equation with 
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Dirichlet type boundary condition^ flT], 

A/ + A/ = Vx G 1], /(x) = 0, Vx G dn, 

where denotes the boundary of the object and A G M is a scalar. There are 
infinitely many pairs {(An,/n)}^i satisfying this equation and they form the 
set of eigenvalues and eigenfunctions respectively. The ordered set of eigenvalues 
is a positive diverging sequence such that < Ai < A2 < • • • < A^ < . . . . This 
infinite sequence is called the Dirichlet spectrum of A^, which we refer simply 
as the "spectrum". In addition, each component of the spectrum is called a 
"mode" , e.g. A^ is the called n^^ mode of the spectrum 

The spectrum contains information on the intrinsic geometry of objects. 
Weyl in 41 showed the first spectrum-geometry link by proving that the asymp- 
totic behavior of the eigenvalues is given as 



n 



2/d 



where Vq is the volume of ft and is the volume of the unit ball in R^. Later 



works, as 27,31,32,37 , extended this result by studying the properties of the 
Green's function of the Laplace operator, and showed that a more accurate 
spectrum-geometry link is given by the heat-trace, which in is given as 

00 00 

Z{t) 4 J2e-'-' = J2a^/2t-'^'+^/\ t>0. (1) 

n=l s=0 

The coefficients of the polynomial expansion, are the components carrying 
the geometric information. These coefficients are given as sums of volume and 
boundary integrals of some local invariants of the shape, 27|37|4Q . For instance, 
as given in 



27 , the first three coefficients are: 



^1/2 = "4(4^)^/2-1/2^^' 
ai = ; -TT777 / K:ddQ, 



where Sq is the surface area (circumference in 2D) and is the mean (geodesic) 
curvature on the boundary of Q. The functional relationship between the eigen- 
value sequence and the coefficients as/2 can be seen in Equation (IT]). This 
connection relates the spectrum to the intrinsic geometry, which is the reason 
why Laplace spectrum is important for the computational study of shapes. 

In addition to the spectrum-geometry link, the eigenvalues of the Laplace 
operator have two other properties which make them useful for shape analy- 
sis, flP. These are: 1) the Laplace operator is invariant to isometric transfor- 
mations and 2) the spectrum depends continuously on the deformations applied 
to the boundary of the object. The advantage of the first property is obvious 



-•^ Other boundary conditions yi eld different eigensystems. Here we are only interested in 
the Dirichlet type. Please refer to pjl for the other types. 
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since isometric transformations do not alter the shape. In addition to this, the 
second property states that there is a continuous link between the differences 
in eigenvalues and the difference in shape, which makes eigenvalues ideal for 
measuring shape differences. 

Unfortunately, it has also been shown that there exists isospectral non- 
congruent objects, i.e. objects with different shape but the same spectrum 15 



Therefore, theoretically the Laplace spectrum does not uniquely identify shapes. 
However, as stated in [3^, practically this does not cause a problem mostly be- 
cause the constructed isospectral non-congruent objects in 2D and 3D are rather 
extreme examples with nonsmooth boundaries. 

The spectral signature, shape-DNA, proposed in [33] is inspired from the 
properties given above. For a given shape 1^, its shape-DNA is the first N 
modes of the spectrum of the Laplace operator defined on Q: [Ai, A2, . . . , Aat]. 
In addition to the properties the shape-DNA inherits from the eigenvalues, the 
authors also proposed several normalisations to obtain almost scale invarianc^ 
The normalisations used in the experiments in 25 , 29 , 33 , 34 are given as A^ 
XnV^^"^ and An An/Ai. 

In 33 , the authors also defined a shape distance based on shape-DNA. 
Either using the original or its scale invariant version, this distance is given as 



N 



1/2 



(2) 



where denotes the object with the spectrum {^n}^i- Using p^^{Qx^Q^)^ 
the authors were able to distinguish between distinct shapes 25 , construct 
shape manifolds based on the pairwise distances and perform statistical com- 
parisons 29 34 . 

However, as also pointed out in 28 , due to the diverging nature of the 
spectrum, p^^ suffers from three essential drawbacks limiting its usability: i) 
differences at higher modes of the spectrum have higher impacts on the final 
distance value even though they are not necessarily more informative about the 
intrinsic geometry, ii) the distance is extremely sensitive to the signature size 
A/", while the choice of this parameter is arbitrary, and iii) the distance cannot 
be defined over the entire spectrum because it does not yield a proper metric 
in that case. Therefore, defining a sound and intuitive distance based on the 
spectrum is still an open question for which we propose a solution in the next 
section. 



3 Weighted Spectral Distance - WESD 

This section presents the proposed spectral distance, WESD, the analysis of the 
heat-trace leading to its definition and its theoretical properties. The structure 
of presentation aims to separate the definition of the distance, which is essen- 
tial for its practical implementation, from the details related to its derivation 
and theoretical properties. In this light, we first present the definitions and 

^We use the term "almost" because scale invariance is an application dependent concept 
and the definition of scale difference between arbitrary objects is a mathematically vague 
notion. A further discussion of scale invariance is outside the scope of this article and we refer 
the reader to [33] , 
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mention the associated properties with appropriate references to the following 
subsections, which contain further details. 

We define the Weighted Spectral Distance - WESD - for two closed bounded 
domains with piecewise smooth boundaries, I^a,^^ C as 



E 



(3) 



with p G M and p > d/2. Unlike the distance given in Equation ([2|, WESD 
is defined over the entire eigenvalue sequence and the factor p is not fixed to 
2. In addition, the difference at each mode contributes to the overall distance 
proportional to lA^ — £,n\/^n^n instead of \Xn — ^n|- The additional A^^^ factor 
(seeming like a simple addition to Equation [2| actually arises from analysing 
the relation between the n*^ mode of the spectrum and the heat-trace, which 
will be presented in Section 3.1 This analysis also provides WESD with a 



geometric intuition. Furthermore, for p > d/2 the infinite sum in the definition 
is guaranteed to converge to a finite value for any pair of shapes. Hence, WESD 
exists. In addition to its existence, WESD also satisfies the triangular inequality 



making it a pseudometric. These points are proven in Section 3.2 Moreover, the 
pseudometric WESD has a multi-scale aspect with respect to p. In Section |3.3 
we show that adjusting p controls the sensitivity of WESD with respect shape 
differences at finer scales, i.e. with respect to geometric differences at local 
level such as thin protrusions or small bumps. Thus, for higher values of p the 
distance becomes less sensitive to finer scale differences. 

In addition to WESD, we define the normalised score for shape dissimilarity 
nWESDd^s 



G [0,1) 



(4) 



W(^]a,^^) 

which maps p{^x^ to the [0, 1) interval using the shape-dependent normali 



sation factor 



C^K- 



2p-i 

d 



The factors C and K are the shape based coefficients defined in Corollary [TJ 
and ("(•) is the Riemann zeta function 42 . Being confined to [0,1), nWESD 
allows us to i) compare dissimilarities of different pairs of shapes and ii) easily 
use the shape dissimilarity in combination with scores quantifying other type 
of differences between objects such as volume overlap in case of matching or 
Jacard's index in case of accuracy assessment. 

One important issue in defining a distance or a score using the entire eigen- 
value sequence is computational limits. In practice we can only compute a finite 
number of eigenvalues and therefore, can only approximate such distances. Con- 
sidering this, here we define the finite approximations of WESD and nWESD 
using the smallest N eigenvalues as 



■ N 

E 

.n=l 
N 



I An — in\^ 
An^n 



1 1/P 



W(nA,f^«) 



(5) 
(6) 
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where is a truncation parameter. Previous works, such as 19,21,22,33 
[35] , also define distances based on finite number of modes. However, their 
view on the distance definition was first to construct finite shape signatures 
and then to define a distance on the signatures. Therefore, the signature size 
was a critical component of the definition itself. Furthermore, the effects of 
the choice of the signature size on the distance values have not been carefully 
analysed in these works. The view presented here defines the distance directly 
using the entire sequence without constructing a finite signature. This alleviates 
the importance of the signature size on the distance. The finite computation 
given in Equations [5] and [6] are viewed as approximations to the distance and 
N as the truncation parameter. In this conceptually different setting, unlike 
previous works, we provide in Section |3.4| a careful analysis of the choice of 
N on the spectral distance. Specifically, we prove that limAr^oo \p{^x^^^) — 
{ftx^^c)\ = and limAr^oo \p{^x^^^) — P^(^A,^^)| = 0. Furthermore, we 
provide a theoretical upper bound for these errors that shows how fast they 
decrease in the worst case leading to a principled strategy for choosing N. 



Section 3.5 ends the section by focusing on the invariance of WESD and 
nWESD to global scale (relative size) differences between objects. Specifically, 
we discuss how an "approximate" scale invariance can be attained for WESD 



and nWESD by following the same strategy proposed in 33 



3.1 Analysis of the Heat-Trace and Derivation of WESD 

We derive WESD by analysing the mathematical link between the spectrum 
of an object and its geometry. This link is given by the heat-trace defined 
in Equation ([T]). Let us consider the heat-trace as a function of both t and 
the spectrum, Ai, A2, ...)• The main question we answer is how much the 

•) function changes when we change the n^^ mode of the spectrum from 
An to ^n- Considering the polynomial expansion equivalent to Z{t^-) given in 
Equation ([T]), one can see that the change in the value •) is directly related 
to the changes in the coefficients as/2 and so to changes in the integrals over the 
local invariants. By analysing the influence of the change in the n*^ mode on 

•), we actually analyse the influence of this change on the integrals over local 
geometric invariants. Following this line of thought, we quantify the influence 
of the change from to on Z{t^-) in terms of and ^n- This can be done 
by deflning 

= / . . . , An-l, An, An+l, . . . ) 

Jo 

— . . . , An-1, ^n, An+1, • • • ) I 

which is simply the Li-norm of the difference between the functions that is 
linked to the difference between An and ^n- Replacing Z{t^-) with its deflnition 
leads to 

POO 

= / \e-^-' -e-^-'\dt (7) 
Jo 

Without loss of generality let us assume ^n ^ ^n- Then 

e-^-^ > e-^-^ for t > 0. 
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We can then evaluate the integral in Equation ^ as 

I — S,n\ 



-Xjit 



captures the influence of the difference at the n*^ mode on Z{t^'). Now, 
aggregating these influences across all modes leads to the definition of WESD 



i/p 



E 



1 1/P 



Surprisingly, the formulation of WESD, which results from the analysis pre- 
sented above, also has very beneficial properties that makes it theoretically 
sound and useful in practical applications. These properties will be analysed in 
the following. 

Before delving into this analysis though let us make two remarks. The first 
relates Psd{'j-) (Equation to the analysis of the heat-trace presented above. 



Remark 1. Let us 



I 

Jo 



oo 



Z (t, An_i, An, An+i, ...) 

Z (t, An-l, ^n, An+1, ...) dt 



Then A^' = A^. Evaluating this integral yields A^'^ 
setting m = 2 psD{-r) can he derived as follows 



oo 



1/2 



This derivation not only relates WESD to psoi'^) ^'^^ ^^^o provides the link 
between psoi'^) the heat-trace. 

The second remark notes the link between WESD and Global Point Signa- 
tures (GPS), a local shape descriptor, presented in 35 . 



Remark 2. GPS, as presented in 13^, is defined for each point in an object 

r 1/2 ^ oo 

Qx as the infinite series GPSnA^) ^ {$A,n(x)} = A-^'7„(x) , where 

I J n=l 

X G cii^d fn{'^) is the n^^ eigenf unction. GPS has a connection to WESD 
arising from the following element-wise integrals 



^A,n(x)^ix: 



a: 



-1/2 



/n(x) 



— A ""^ , 



where the equality arises from the fact that eigenf unctions form an orthonormal 
basis in Vtx ll^ i.e. J^^ fn{'x.)fm{^)dx = 5{n — m) with 5{-) being the Dirac^s 
delta. Consmering this integral, WESD can also be regarded as a distance be- 
tween GPS^ of two objects as 



E 



*A,n(x)dx- / $|,„(x)cZx 
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This link also provides an alternative view on the normalisation factor An 
used in GPS. In [35j author justifies this normalisation factor by noting that for 
an object the Green's function can be written as an inner product in the GPS 
domain, see Section 4 in fS^ - This is later used to argue the geometric meaning 
of GPS as authors point out the use of Greenes function in different shape pro- 
cessing tasks. Our link between GPS and WESD provides an alternative view 
on the normalisation factor as it connects this local signature to the heat-trace 
Z{t). 



3.2 Existence of the Pseudometric WESD 



WESD is defined as the limit of an infinite series as given in Equation (|3|. For 
such a distance to be a proper one, actuahy a pseudometric in this case, the hmit 
of the infinite series should exist for any two spectra. In the case of WESD, this 
is not evident because it is defined over the entire spectra and each spectrum 
is a divergent sequence. The first corollary presented below proves that when 
p > d/2 WESD indeed satisfies this condition, i.e. the infinite series converges. 
The corollary further provides an upper bound for this limit, which is used to 
construct nWESD. We would like to note that for the ease of presentation, the 
proofs for all the following corollaries and lemmas are given in Appendix B in 
the supplemental material. 



Corollary 1. Let Qx C and 

piecewise smooth boundaries and {A}^ 
Then the weighted spectral distance 



C be any two closed domains with 
i and {^}^i be their Laplace spectrum. 



oo 

E 

n=l 



I An — 



1/p 



converges for p > f . Furthermore, 



p{nx,n^) <{c + K- 



(8) 



where ("(•) is the Riemann zeta function and the coefficients C and K are given 



as 



C 



K 
V 



E 

i=l,2 

d^ 



BdV 



1 



1 

d 



i-l 



_d-47r2 V ; II (i + 2.64_ 
max(V(l]A), II = max(Ai,^i), 

where V{-) denotes the volume (or area in 2D) of an object. 



The Inequality (|8| states that WESD has a shape-dependent upper bound. 
We thus can map the WESD to the [0, 1) interval through normalising it with 
this upper bound. The nWESD score, given in Equation [4] is constructed based 
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on this strategy. Since its existence is established next we prove that WESD 
is a pseudometric, i.e. satisfies the other criteria to be a pseudometric, such as 
the triangle inequality. 

Corollary 2. p{^x, is a pseudometric for d>2. 

We note that WESD is not a metric because the spectrum is invariant to 
isometrics, which is a desirable property for shape analysis. However, in addition 
to this, the spectrum is also invariant to isospectral non-congruent shapes. This 
is not desirable but does not cause problems in practice as discussed in Section 
[2] and also confirmed in our experiments. 

3.3 On the multi-scale aspect of WESD 

The previous section highlighted the role of p on the convergence properties 
of WESD and therefore on its existence. We now demonstrate that p also 
provides WESD a multi-scale characteristic. The sensitivity of WESD to the 
shape differences at finer scales depends on the value of p. Specifically, we show 
that the higher the value p the less sensitive WESD is to finer scale details and 
its sensitivity increases as p gets lower. 

The multi-scale aspect of WESD arises from the relationship between the 
Laplace operators and heat diffusion processes 12 . We first present an intuitive 
summary of this relationship, which is about the multi-scale aspect of Z{t) and 
t in particular. For a more mathematical treatment we refer the reader to 38 . 
As stated in 38 and 28 , t can be interpreted as the time variable in a heat 
diffusion process within an object. A useful visual analogy to consider here is 
the Laplacian smoothing of a surface where t would correspond to the amount of 
smoothing. Similar to the surface smoothing, as t increases, the local geometric 
details of an object, such as sharp ridges or steep valleys, lose further their 
infiuence on the Z{t) value. As a result Z{t) becomes somewhat insensitive to 
these local geometric details, in other words shape details at finer scales. From 
an alternative view, the value of Z{t) loses its information content with regards 
to local geometric details. This effect intuitively summarizes the multi-scale 
characteristic of the heat-trace with respect to t. 

Having explained the multi-scale aspect of Z{t)^ we now analyse how this 
aspect is refiected upon the eigenvalues. To do so let us define the influence 
ratio V{n^t) = ^z{t) ' This ratio captures the infiuence of the n*^ mode on 
the heat-trace. In other words, the higher the ratio, the higher the infiuence of 
\n on the value of Z(t) at that specific t. The following lemma compares the 
influence ratios of different modes and how this comparison depends on t. 

Lemma 1. Let Q.x C represent an object with piecewise smooth boundary 
andV(l^t) = ^z{t) corresponding influence ratio of mode I att. Then for 

any two spectral indices m > n > 

V{n,t) >V{m,t), yt>0 

and particularly for two t values such that ti > t2 

V{mM) V{m,t2) 
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The first inequality of tiie lemma indicates that the lower modes in the 
spectrum have more influence on the value of Z(t) than the higher modes. The 
second inequality shows that the influence of the higher modes become more 
prominent as t decreases. Considering that for lower t values Z{t) is more 
informative with regards to shape details at finer scales, Lemma [T] suggests that 
the higher modes are more important for finer scales than for coarser scales. We 
illustrate this observation on a synthetic example shown in Figure [2] with the 
pair (a) + (b) being an example showing coarser shape differences and the pair 
(a) + (c) showing finer differences. The plots given in Figure [2] (d) and (e) show 
the corresponding spectral differences observed at modes between 1 and 150. 
Between (a) and (b) the shape differences are at the coarse level. According to 
Lemma [1] these differences should show up at the very first modes. On the other 
hand, between (a) and (c) the differences are at a finer scale and furthermore 
the objects are very similar at the coarse level. Lemma [T] states that these 
differences therefore, should show up at higher modes and the differences at the 
lower modes should be low. Satisfying these expectations, the differences at the 
first few modes shown in plot (d) have relatively large values compared to the 
ones in plot (e). Furthermore, the amplitude of the differences at higher modes 
are generally larger in plot (e) than in plot(d), especially after 100. 

In order now to connect these findings to WESD and p let us present the 
following corollary, which studies the influence of p on the components inside 
the infinite sum defining the distance. 

Corollary 3. Let Qx and be two objects with piecewise smooth boundaries. 
Then for any two scalars with with p > d/2, q > d/2, p > q and for all n with 
\^n — in\ > there exists a M > n so that Vm > M 




Thus, the relative contributions of the higher spectral modes on p{Vtx.,Vt^) 
with respect to the contributions of the lower modes depend on the value of p. 
Specifically, the higher spectral modes become more influential as p decreases. 
Combining this finding with the result of Lemma [l] we follow that as p increases 
WESD gives less importance to differences at higher spectral modes and there- 
fore becomes less sensitive to the shape differences at finer scales. This provides 
WESD with a multi-scale aspect with respect to p and also provides us the 
intuition for choosing a proper value for p. 

3.4 Finite Approximations of WESD and nWESD 

One of the important practical questions regarding spectral distances is the num- 
ber of modes to be included in the calculation of the distance. The computation 
of eigenvalues and eigenfunctions can be expensive and inaccurate especially for 
the higher modes. Therefore, spectral distances require the user to set a finite 
number of modes to be used. This parameter is often referred to as the signa- 
ture size. Having defined the distance over the entire sequence, we refer to it 
as the truncation parameter. This actually provides a different perspective on 
the number of modes used to compute the distance. In previous works, such 
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50 100 150 

modes 

(d) between (a) and (b) 




50 100 150 

modes 

(e) between (a) and (c) 



Figure 2: Multi-scale characteristics of different spectral modes: (a), (b) and (c) 
show three synthetic shapes. In (d) we plot the absolute differences between the 
corresponding modes of (a) and (b) with respect to the spectral index. In (e) 
we plot the same difference for the shapes in (a) and (c). The shape difference 
between (a) and (b), which is at a coarser level, is already captured at the lower 
spectral modes. The difference between (a) and (c) results in lower differences 
in lower spectral modes because these objects are more similar at a coarser level. 
At the higher spectral modes though the difference between (a) and (c) becomes 
more prominent since these objects differ more substantially at the finer scales. 
The plots in (d) and (e) demonstrate that the higher modes for a given object 
are more important for finer scale shape details. 



as 21,33,35 , the value of this parameter, viewed as the signature size, is often 
set arbitrarily and its effect on the distances have not been carefully analysed. 
Here, viewing it as a truncation parameter, we study its influence. Specifically, 
we formulate the difference between using the entire spectra to only using a 
finite number of modes as an approximation/truncation error. So we analyse 
how this error changes with respect to the truncation parameter. We specifi- 
cally show in the next corollary that the errors in approximating WESD and 
nWESD by the first N modes converges to zero as N increases. Furthermore, 
we provide an upper bound for both errors as a function of N. 

Corollary 4. Let p^{Qx^ Q^) be the truncated approximation of p{flx, fl^) based 
on the first N modes and p^{Qx, Q^) ofp{Qx^ ft^). Then \/p > d/2 

lim \p-p^\=0 
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N 



lim \p 



0. 



Furthermore, for a given N >3 the truncation errors \p — and \p — can 
be bounded by 



\p-p^\ < IC^K 




(9) 



\p-Pn\ < 1 



[Zln=3 in) 



C^K- 



(10) 



The above corollary has important practical implications. First of all, the 
sensitivities of p^ and p^ with respect to TV decreases as TV increases. For any 
application relying on the shape distances, such as constructing low dimensional 
embeddings, this reduced sensitivity is particularly important as it provides 
stability with respect to N both for the distance and for the application using 
the distance. We note that the opposite is true for p^^, which is one of the 
main disadvantages of this distance. 

In addition. Corollary [4] can guide the choice for the number of modes 
and the norm type p. Specifically, the error upper bounds given in Equations |9] 
and [To] provide the worst case errors for a given pair of shapes without the 
need to compute the eigenvalues. So for instance, once a number of modes are 
computed then based on the distance value obtained so far and the worst case 
error computed using the upper bounds, one can decide whether to compute 
more modes or not. Furthermore, these upper bounds are shape-specific as 
they depend on C and K. One can go one step further and define a shape- 
independent residual ratio for N > 3 and p > (i/2 as 



c(f)-i-G) 



(11) 



that satisfies R{N,p) > p — p^ for which the proof is given in Proposition 1 
in Appendix B. Based on this, R{N^p) can be used to select the parameters N 
and p as it quantifies the quality of the approximation for a given pair of (iV,p) 
in terms of the error upper bounds. 

In Figure [sj we plot R{N^p) versus N for different settings of p and d = 2, 3. 
Besides the obvious point that the error upper bound decreases for increasing 
N we also notice that i) the behavior in 2D and 3D are similar and ii) the rate 
of decrease of the error upper bound is much faster for higher p. Considering 
the multi-scale aspect of WESD captured in p, this behavior is interesting. It 
demonstrates that the choice of p and N are correlated and suggests a trade-off 
between the rate of decrease of the truncation error and the sensitivity of WESD 
to shape differences at finer scales. In theory, the choice of these parameters 
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Figure 3: Choosing N: The figures plot the residual ratio R{N,p) versus N for 
different p values in 2D (left) and in 3D (right). As expected the error upper 
bound drops with increasing N. The rate of decrease also becomes faster with 
increasing p. This inverse relation suggests the trade-off between N and the 
sensitivity of WESD to finer scale shape differences since WESD becomes less 
sensitive as p increases, see Section 13.3 



depends on the application and the expected shape differences. If one expects 
coarse scale differences then choosing a large p and small N might be sufficient. 
However, if one is interested in finer scale differences then a small p value will 
be required, which in turn will require a large N value to have a decent approx- 
imation. The important aspect of R{N^p) is that it is universal, i.e. it does not 
depend on the objects. So it can be used in any type of application to choose the 
parameter pair p and to have a rough estimate of the computational costs for 
computing the distance WESD. We note once again, the specific values should 
be chosen based on the application and the shapes at hand. 



3.5 Invariance to global scale differences 

We end this section studying the impact of global scale differences on WESD 
and how invariance to such differences can be attained. We would like to note 
that the notion of global scale in this section refers to the relative size of an 
object, which is not to be confused with the notion of multi-scale used in Sec- 
tion l33| The spectrum of an object depends on the object's size, i.e. a global 



scale change alters all the eigenvalues by a constant multiplicative factor 11 
As a result, the global scale difference between two objects contributes to the 
spectral shape distance WESD. In some applications this contribution might 
not be desirable, for instance in an object recognition task, where objects in the 
same category have varying sizes. Therefore, it is a useful property of a shape 
distance to allow invariance to global scale differences. 

Renter et al. [33^ proposed different approximations for normalising the 
effects of scale differences on the spectrum. In particular, the authors use two 



different normalisations in their experiments in [25l[29l[34 . Both normalisations 



directly act on the eigenvalues. The first one normalises the eigenvalues with 
respect to the volume (area in 2D or surface area for Riemannian manifolds) 
and is given as ^uVq^^' The second one normalises the eigenvalue with 

respect to the first eigenvalue in the sequence as A^/Ai. Both of these 

strategies can be used when computing distances with WESD. Furthermore, 
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since these strategies do not alter the mathematical characteristics of the entire 
spectrum the theoretical properties of WESD and nWESD hold either way. For 
our experiments we adopt the first strategy, volume normalisation, using the 
volume as defined in Euclidean geometry. When using the volume normalised 
eigenvalues, the only change that applies to the technical details presented so 
far is V in Equation [s] becomes V = 1. The rest applies directly without any 
modification. 

We would also like to note that estimating the global scale difference between 
two arbitrary objects is not always a well-posed problem. It is especially hard 
when the objects are of different category, e.g. an octopus and a submarine. 
Furthermore, the scale normalisation is application dependent and it might not 
be desirable for all applications. In Section 4.3.2| we present such an example 
where we analyse the temporal change of the left ventricle shape during a heart 
cycle. In this case, the volume change is essential for analysing the heart of the 
same patient so that scale invar iance is not appropriate. 



4 Experiments 



This section presents a variety of experiments on synthetic and real data high- 
lighting the strengths and weaknesses of WESD and nWESD. We start by briefiy 
explaining the details of the numerical implementation of WESD used in the 
experiments presented here. Then in Section |4.2[ the proposed distances are 
applied to synthetically generated objects demonstrating that 

(i) Ordering objects with respect to their shapes using nWESD results in a 
visually coherent series (Section 4.2.1), 



(ii) WESD is useful for constructing low dimensional embeddings, in partic- 
ular it yields stable embeddings with respect to the signature size 



(Section 4.2.2) and 



(iii) WESD is a suitable distance for shape retrieval, which is shown through 
experiments on the SHREC dataset 25 (Section 4.2.3). 



Lastly, in Section [43] WESD is applied to real objects extracted from 3D medical 
images. We focus on two examples from a wide variety of applications WESD 
and nWESD can be beneficial to: population studies of brain structures and 
analysis of 4D cardiac images. 



4.1 Implementation Details 

There are two different aspects in the implementation of WESD: the numeri- 
cal computation of the Laplace spectra and the parameter settings. First, any 
numerical method tailored towards computing the eigenvalues of the Laplace 
operator can be used to compute WESD. Examples of such method are listed 
in [l][33]. Our specific implementation represents objects simply as binary im- 
ages with the foreground defining ft. Using the Cartesian grid of the image, it 
discretizes through finite difference scheme (see also Chapter 2 of fT). This 
step yields a sparse matrix of which we compute the eigenvalues via Arnoldi's 
method presented in |2j and implemented in MATLAB®. We choose this spe- 
cific implementation as 1) it is simple 2) it does not introduce any additional 
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parameters and 3) when working with images it avoids any extra preprocessing 
steps, such as surface extraction or mesh construction. 

With regards to the second implementation aspect, we set the parameters N 
and p empiricahy. Based on Section [3^ we set p = 1.5 in 2D and p = 2 in 3D. 
These values result in a relatively fast diminishing upper bound of the truncation 
error with respect to N (see Figure |3| while being sensitive to shape differences 
at finer scales. In both 2D and 3D, we chose N = 200 for the number of modes 
as the truncation error seemed to vanish at that point. Furthermore, in addition 
to the theoretical considerations on the effects of and p on WESD given in 
Section 3.4, in Sections 4.2.2 and 4.2.3 we experimentally study the effects of 
these parameters on applications using WESD, specifically on constructing low 
dimensional embeddings and shape retrieval. 



4.2 Synthetic Data 

We conduct three experiments: first two are on 2D objects and the last one is 
on 3D objects. For all of the experiments, we use the scale invariant versions of 
the spectra obtained by normalising the eigenvalues with the object's volume as 
described in Section 13.51 As a result the distances WESD and nWESD become 
"almost" invariant to global scale differences. 



4.2.1 Ordering of Shapes 

For the first experiment we created two synthetic datasets. Each dataset consists 
of a reference object and random deformations of this reference. These deformed 
versions are generated by transforming the reference via random deformations 
of varying magnitude and amount of nonlinearity. As a result the datasets 
contain objects that are very similar to the reference ones and objects that are 
substantially different. Figures |4ja) and (b) show some examples from these 
datasets where the binary images to the very left show the reference objects. In 
the first dataset, the reference object is a disc and in total there are 500 random 
deformations of this reference disc. The first 400 are generated via non-linear 
deformations while the last 100 are isometric transformations. In the second 
dataset, the reference is a slightly more complicated object (see Figure |4jb)) 
and in total there are 400 random transformations of this reference. The first 
300 are generated by non-linear deformations and the last 100 produced via 
isometric transformations. All objects are discretized as binary maps with a 
size of 200 x 200 pixels. The numerical computations are performed on these 
image grids as discussed earlier. 

We computed the nWESD scores {p^ with p = 1.5 and N — 200) between 
the reference and the deformed objects in each dataset. Based on these scores, 
we then ordered the deformed objects according to their similarity in shape to 
the reference. Figures |4] (c) and (d) show examples of the resulting orderings. 
We notice that the orderings are visually meaningful , i.e. the further the 
deformed objects visually deviate from the references, the higher their nWESD 
score is. Furthermore, all the objects generated via isometric transformations 
yielded scores close to zero as a result of the invariance of the proposed scores 
to this type of transformation. 
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Figure 4: Shape-based ordering of objects: We generate two artificial datasets 
each consisting of a reference object and its random deformations. Samples 
from the datasets are shown in (a) and (b). The binary images to the very left 
show the reference objects for each dataset. We then ordered all the deformed 
objects with respect to the nWESD scores between the object and the reference. 
The graphs in (c) and (d) plot these orderings. Based on visual inspection the 
ordering is quite reasonable. 



4.2.2 Low Dimensional Embeddings 

In the second experiment we focus on creating low dimensional embeddings. 
We compare the embeddings constructed by WESD with the ones constructed 
using (Equation ([2|), the distance proposed in 
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We do so based on 



the TOSCA dataset (toolbox for surface comparison and analysis), [6j[7]. This 
dataset contains binary segmentations of 5 human, 5 centaurs and 5 horses 
as shown in Figure [sja). We compute the pairwise affinity matrices between 
objects via pf^ and WESD {p^ with p = 1.5). We then apply the ISOMAP 
algorithm 39 to these matrices, which maps the 15 objects to a 2D plane based 
on the pairwise shape distances. We repeat this experiment for affinity matrices 
computed using different number of spectral modes, i.e. N = 50, 100, 200, to 
demonstrate the effect of the signature size (truncation parameter) on both 
distances. 
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Figure 5: Low dimensional embeddings: (a) The 15 objects used in this experi- 
ment. The graphs plot the 2D embeddings of the objects based on the affinity 
matrices constructed by p^^ and WESD (p^). Each row presents the results 
based on different number of eigenvalues: 50, 100 and 200 from top to bot- 
tom respectively. The structures of the 2D embedding based on p^jj are quite 
different for different N. WESD however, produces embeddings that are simi- 
lar. This demonstrates the stability of the embedding with respect to when 
WESD is used. 
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The plots in Figures [5|b),(d) and (f) present the resulting 2D embeddings of 
the dataset using p^^. The embeddings are substantially different for different 
N. This variation arises due to high sensitivity of p^^ towards the signature 
size N. Moreover, the embeddings obtained using higher N are less satisfac- 
tory in terms of separating the three different object classes. This is actually 
as expected since the spectral modes with higher indices dominate the value 
of pg]j even though they are not informative with regards to the overall ge- 
ometry and thus, negatively impact the outcome. The plots in Figures [sjc), 
(e) and (g) present the embeddings obtained using WESD. Apart from simple 
coordinate flips (arising from ISOMAP's indifference to signs) the embeddings 
obtained at different N are very similar. This shows that the construction of 
the low dimensional embedding is stable with respect to N when WESD is used. 
This is a direct consequence of the convergent behavior of WESD discussed in 



Sections 3.2 and 3.4 As illustrated by this experiment, this property has very 



important practical implications. 



4.2.3 Shape-based Retrieval of 3D Objects 

In this last experiment with synthetic data, we focus on the application of shape- 
based object retrieval^ i.e. given a test object identifying other "similar" objects 
within a dataset using shape information. Similarity in this context can be de- 
fined in various ways but the definition used here is semantic similarity, meaning 
that objects that are of the same semantic category (e.g. human bodies, aero- 
planes, etc) are similar and objects of different categories are not. Shapes of 
similar objects have similar traits and properties. Shape distances used for re- 
trieval purposes should be able to capture these traits yielding the lowest values 
between similar object pairs. Here, WESD's value for shape-based retrieval is 
evaluated using the publicly available dataset SHREC presented in 25 . 

SHREC dataset consists of 600 3D non-rigid objects from 30 different cat- 
egories, i.e. 20 objects per category. Objects from the same category differ 
with substantial non-linear deformations, which makes retrieval in this dataset 
challenging. To evaluate the retrieval accuracy of WESD, first each object was 
converted from its original watertight surface mesh discretization to a 3D binary 
image using the Iso2mesh software packag^ Then pairwise shape distances 
across the entire dataset were computed using WESD and the 600 x 600 affinity 
matrix was constructed, where each entry is a pairwise distance. This affin- 
ity matrix was then evaluated using the software provided with the dataset ^ . 
The evaluation consists of a variety of retrieval accuracy scores such as Nearest 
Neighbor (NN), First-Tier (FT), Second-Tier (ST), E-Measure (E), Discounted 
Cumulative Gain (DCG) and Precision-Recall curve. The first two rows of the 
table in Figure [6]^ a) list these scores obtained using WESD for two different set- 
tings of the p and N values. Additionally, the last three rows of the same table 
show the results obtained using p^^ (Equation |2j 33 ), as listed in 25 us- 
ing two different types of scale normalisation (norml: normalising with respect 
to the first eigenvalue, normA: area normalisation, see Section |3.5| for further 
details). These accuracy scores show that WESD and pgj-, perform very simi- 

^Available at http: / / www.it 1. nist .gov / iad / vug / shar p / contest /201 1 /NonRigid / 1 
^http://iso2mesh. sourceforge.net/cgi-bin/index.cgi 

^We note that for these latter results a slightly different notation is used here than in [25] 
to conform to the overall notation of this article. 
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lar in retrieval from the SHREC dataset. Furthermore, Figure [6[b) shows the 
precision-recah curve of WESD {p = 3.15 and = 100) for the entire dataset. 
The curve is very similar to the best curve obtained using pgj^ shown in 
Once again, this confirms that both distances perform similarly. 

Lastly, the graphs shown in Figur^c) and (d) provide an analysis of the 
retrieval results with respect to the parameters N and p. Graphs in Figure [6]^c) 
plot the change of different retrieval scores with respect to the number of modes 
used i.e. signature size, keeping p fixed at 3.15. Graphs in Figure |6jd) plot 
the changes with respect to the norm type p keeping N fixed at 100. These 
graphs show that as N increases the scores seem to increase slowly and then 
converge. On the other hand, p has a stronger effect on the results than 
particularly on FT, ST and E scores. However, the changes in the scores with 
respect to changes in TV or p are rather small especially compared to the rela- 
tively larger fluctuation of the FT score of p^^ with respect to the two sample 
values provided in the table in Figure [6]^a). 

The experiment presented above showed that the retrieval power of WESD 
is similar to that of the distance p^jj proposed by Renter et al. |33 . The 
soundness and theoretical properties of WESD do not come at the expense of 
lower retrieval power. On the contrary, WESD is able to leverage the descriptive 
power of the spectra while its properties guarantee that it does not suffer from 
similar drawbacks as other distances, such as sensitivity to signature size. 



4.3 Real Data 

The experiments on real data are conducted on segmentations of 3D structures 
obtained from magnetic resonance images (MRI). First, we apply WESD to 
subcortical brain structures. The experiment demonstrates WESD's capabilities 
to differentiate categories of objects even in the presence of high intra-class 
variability. In the second experiment, we focus on temporal analysis of cardiac 
images. We apply nWESD to delineations of the blood pool of the left ventricle 
obtained from 3D + time cardiac MRI. The experiment shows that the shape 
dissimilarity measurements between time points correlates with the dynamic 
processes of the beating heart. 



4.3.1 Clustering Sub-Cortical Structures 

Medical research frequently relies on morphometric studies analysing anatom- 
ical shapes from medical images [2]. In this experiment we construct a low 
dimensional embedding of subcortical structures extracted from Magnetic Res- 
onance Image (MRI) scans of different individuals based on WESD as well as 
shape-DNA based distance, p^^^ as proposed in [ss]. 

For this experiment, we use the publicly available LPBA40 dataset [36 Q The 
dataset contains manual segmentations of various subcortical structures from 
MRI brain scans of 40 healthy subjects. Figures [t]^ a), (b) and (c) show some 
examples from these structures. The are two main difficulties associated with 
such datasets. First, the structures have very large intra-class (inter-subject) 
variability, i.e. the shape of an anatomical structure is often very different across 
subjects. Second, the segmentations were obtained by manually delineating 
the 3D objects on successive 2D slices. This creates inconsistencies between 

^website: http: / / www. loni.ucla.edu / Atlases /LPB A40 
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segmentations in two successive slices. Such inconsistencies in the end manifest 
themselves as local artefacts on the object. The protrusion that can be seen 
on the top of the second hippocampus in Figure [7|^a) is an example of such an 
artefact. These artefacts can influence shape distances negatively. 

We select six structures for each patient: left /right caudate nucleus, left /right 
putamen and left/right hippocampus, resulting in 240 structures in total. We 
then create pairwise affinity matrices of the 240 structures ffist using pgj-, with 
N = 200, as proposed in [33^, and then WESD (p^ with p = 2 and N = 200). 
Finally, we use the ISOMAP algorithm [39] to construct 2D embeddings of the 
structures. Figures and (e) show the resulting embeddings. We observe 
that the embedding obtained via WESD well clusters the data with respect 
to the anatomical structures. The separation of the clusters for the SD case, 
however, is more ambiguous, especially between putamen and hippocampus. 

The embeddings presented above were obtained by directly using the manual 
segmentations without any preprocessing. A natural question is how do these 
embeddings change if the effects of various artefacts are reduced say via surface 
smoothing. To answer this question, we smooth the surface of the anatomical 
3D models and recomputed the embeddings, which are shown in Figures 
and (g). The embedding obtained with p^^, although to a lesser extent, still 
suffers from similar ambiguity as in Figure ^d) . The new embedding based on 
WESD on the other hand, compared to Figure [Tf^e), even more clearly sepa- 
rates different anatomical structures. However, we also note that this type of 
preprocessing can also produce undesirable artefacts such as altering the topol- 
ogy of the anatomical object. This is the case for one caudate in Figures 
and (g), which ends up as an outlier that is clearly separated from the other 
data points. Considering this, the fact that WESD is able to produce visually 
pleasing embeddings (see Figure [7|^e)) without the need of preprocessing is an 
advantage. 

4.3.2 Analysing Heart Function in 4D MRI 

Four-dimensional imaging of patient anatomy is gaining interest in the medical 
community. The temporal analysis of anatomical structures is used to extract 
the characteristics of related dynamic processes, which often indicate certain 
pathologies 3,16,26 . In this section, we apply nWESD to the shapes of the 
hearts extracted from four dimensional cardiac images of five different patients. 
The scan of each patient captures a full cycle of one heartbeat as a series of 
20 3D images. Each image shows the left ventricle (LV) at a specific point in 
the cycle, from which we manually segment the corresponding blood pool. Our 
reference is the blood pool extracted from the ffist frame (diastole). We compute 
the nWESD scores between this reference and all other shapes extracted from 
the series of images. Here, we do not normalise the eigenvalues with respect 
to the global scale since size change is an important aspect of the heartbeat 
dynamics. The graph given in Figure |8] shows the results of these measurements 
over time across the five patients. The figure also shows some exemplary images 
and shapes. We observe that the symmetry of the heartbeat along the systolic 
(as the blood pumps out of the LV pool) and the diastolic phases (as the blood 
fills in the pool) is well captured with the nWESD score. Furthermore, the end- 
systolic phase (the time point with the largest distance w.r.t. the reference) 
is at different time points for different patients, which is to be expected since 
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the different patient scans are not synchronized in time. In summary, WESD 
weh captures the dynamics of the beating heart, which is to be expected given 
the continuous hnk between the differences in eigenvalues and the difference in 
shape (see Section [5]). 



5 Conclusion 

This article proposed WESD, a new spectral shape distance defined over the 
eigenvalues of the Laplace operator. WESD is a theoretically sound shape met- 
ric that is derived from the heat-trace. The theoretical analysis given in this 
article presented and proved the properties of WESD related to its existence, 
comput ability and multi-scale aspect. The presented experiments showed that 
the theoretical properties of WESD have many practical advantages over previ- 
ous works. These experiments further highlighted that WESD is beneficial for 
various applications. 

A Bounds on the Laplace spectrum 

Li and Yau in [2? proved that the Laplace spectrum has the following universal 
lower bound 

We notice that this lower bound does not depend on the shape. 



A. 0.3 Upper bounds 



Several authors have investigated the upper bounds and the relative growth 
rate of the eigenvalues of the Dirichlet spectrum 32* .In 43 , Yang provides an 
upper bound for the growth rate of the components for n > 1 as 



1 

n ^ — ^ 



(13) 



This equation can be transformed into a sequence of upper bounds by only 
knowing the first eigenvalue Ai. Although sharp for the first few eigenvalues, 
the upper bound is too relaxed for the remaining modes. Cheng and Yang 
in 10 provides a much sharper upper bound for larger values of n and is valid 
for n > 2. 



where 



An+i < Co{d,n)ndXi, 

a(min((i, n — 1)) 



(14) 



Co{d, n) 



and a{p) = 2.2-41og(l 



d 

a(l) < 2.64 a(2) < 2.27 



50 



-) for p > 3, 



where the bound only depends on the first eigenvalue and furthermore it is 
consistent with Weyl's asymptotic growth law. 
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B Proofs of Lemmas and Corollaries 

Before presenting the proofs for the coroharies let us provide a lemma that will 
be useful throughout this section. 



Lemma 2. For any a, 6 G M such that a > b > the function f{a^ b) 

ah 

increasing b. 

vof. Sin 

^ and ^ 



^ ^ increases monotonously with increasing a and decreases monotonously with 



Proof. Since / is different iable it suffices to look at its partial derivatives ^ = 



db 



1_ 
"62 • 



□ 



Corollary [l| Let C and ft^ C be any two closed domains with 
piecewise smooth boundaries and {X}'^^i and {^}^i be their Laplace spectrum. 
Then the weighted spectral distance 



converges for p > Furthermore, 



E 



I An — I ^ ^ 



i/p 



p{nx,Vt^) <{C + K- 



(15) 



where (■(•) is the Riemann zeta function and the coefficients C and K are given 
as 



C 



E 



d+2 



BdV 



K 4 
V ^ 



=1,2 

d + 2 



d • 47r2 \^ 

(B.vf-i. 



1 / d 
li'\d^^ 

d 



i-l 



d-^TT^ V ; fi (i + 2.64_ 

max{V{nx)^V{n^)), fi = max(Ai,^i), 

where V{-) denotes the volume (or area in 2D) of an object. 



Proof. The following inequality results from combining the bounds specified in 
Section [a] with Lemma [2] 



Xn-U ^ ^ + 2 I BdV 



for n = 1,2 and for n > 3 

I An — £.n\ 



n 



ji \d^4. 



An^n 



< 



< 



d 


+ 2 


(T) 


'-1 


1 


d- 


47r2 




/i 




d 


+ 2 


(T) 


2 


d 


d- 


47r2 






{d^2.QA)n-d 
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Based on this component- wise bound we can write the infinite sum without the 
first two terms as 



E 

n=3 



n sn 



n=3 ^ 



which for p > | converges to 



n=3 



C 



f2p 
V d 



and diverges for p < |. Consequently, /9(riA,^^{) converges for p > |. Further- 
more, extending the sum with the upper bounds forn = 1, 2 the following upper 
bound for the distance between fi^ and $7j holds 



□ 



Corollary [2} piflx^ft^) is a pseudometric for d > 2. 
To ease notation we define 



This leads to 



Proof. Vl^A C M^, C The first three points for this proof are triviah 

- p{rtx,rt^) > since ^^(^a,^^) > Vn. 

- p(l^A, ^a) = since |An — A^l =0 Vn. 

- p{nx, ^a) since \Xn - = - AnI Vn 

In order to prove the triangle inequality let us proceed with the case A^ > ^n- 
The inverse case follows exactly the same way Now, let C be an arbitrary 
closed domain with piecewise smooth boundaries whose spectrum is given as 
{^nl^i- Investigating 

I An -Vn\ ^ \Vn - 

we notice that for each n there are only three possible cases: 
1- An > r/n > ^n, for which 
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2. An > ^ ^n, for which ^^(^A^^ry) ^ Qn{^x^^^) as a result of the 
Lemma|2] Due ^^(^^,^77) > 0: 

3. T]n ^ Xn ^ ^n, foi" which ^^(^r;,^^) ^ ^n(^A7^^) as a result of the 
Lemma [2] once again. And as in the previous case, due to Qni^'q^ ^a) ^ 
we have 

Thus Vn ^n(^A7 ^77) + Qni^T]^ > ^n(^A7 ^^)- Siucc p > 1 as p > | for (i > 2 
the Minkowski Inequahty states 



oo 



> 



E (^n(^A, ^77) + ^n(^r7, ^^)f 



1/p 



When combined with the previous results, the outcome is the triangle inequality: 

□ 

Lemma [l} Le^ ^7^^^^ represent an object with piecewise smooth boundary and 
V(l^t) = ^z{t) corresponding influence ratio with respect to mode I and 

t. Then for any two spectral indices m > n > 

V{n,t) >V{m,t), Vt>0 

and particularly for two t values such that ti > t2 

V{mM) V{mM) 

Proof. The proof follows the properties of the exponential function and the 
properties of the spectrum of the Laplace operator. For n < m we know that 
An < Am which leads to e"'^^^ > e"'^'^* Vt > 0. Since the denominators are the 
same for both P(n, t) and P(m, t) then 

V{n,t) > V{m,t) yt > 0. 

For the second part of the lemma, we first compute the ratio 



V{m,t) 
V{n,t) 



-iXm-Xr,)t 



Now based on Am > An and e ^^"^ '^'^^^ is monotonously decreasing with in- 
creasing t, it follows for ti > t2 that 

^jrn.ti) ^ ^_(A^_A^)ti ^ g-(A^-A.)^2 _ ^(^^^2) 



V{nM) 



V{nM) 



□ 
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Corollary [sj Let Qx and be two objects with piecewise smooth boundaries. 
Then for any two scalars with p > d/2, q > d/2, p > q and for all n with 
\^n — in\ > there exists a M > n so that Vm > M 



w - ( 

Proof From Corollary [l] we know that the series 

l^m I 



E 



converges. Then based on Cauchy's convergence criterion for series 

'lA 



lim 



0. 



In other words, Ve > there exists a M such that 

'lA 



Am ^77 



< e, Vm > M. 



Let n be an arbitrary index such that \Xn — £,n\ > 0- Consequently also for 
I An — ^n|, there exists a M such that Vm > M 

Since p > q we can find a /c > 1 such that p = kq. Then based on the above 
inequality Vm > M 



/ I Arn, — ^TTT. I ^ ^ 

/ |A„-t„l V 



< 



/ \^m—^m I ^ 

/ I Am, — ^TTT, I ^ ^ 



□ 



Corollary 51 Let p^{Qx,Q^) and p^{Qx,fl^) be the truncated approximations 
of p{^x^ (I'^d p{^x^ respectively, using the first N modes. Then 



lim \p-p^\=0 



N 



lim \p — p 



N\ 



0. 
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Furthermore, for a given N >3 the truncation errors \p — and \p — can 
be bounded by 



p-p^\ < {C^K 



C^K 



1 \ ^ 



E 

71 — ?> 



\p-Pn\ < 1-s 



C + K- 



cm-i-G)^ 



Proof. As before, to ease notation, let us again define 

^nsn 



(16) 



Then based on Corollary [T] we know that the sum Yl^=i On exists and thus also 
the partial sums converge 



lim 

N^oo 



N 



n=l n=l 



N 



= Ji^^^E^»-E^» = o- 



n=l n=l 



Based on 



Va, 6, G M with a, 6, d > and a'^ - 6^ ^ ^ a - 6 ^ 
and 

we reach 



lim \p — pn\ = lim 

N^oo N^oo 



N 



As the denominators for both and p are the same, the above limit also yields 
limAT^oo |p- PatI = 0. 

The upper bounds for the truncation errors now is a direct result of Corol- 
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lary [T] as 



\P- P 



N 



< IC + K- 



C + K- 



n=l 



■ AT 

E 



2p 



\P-Pn\ 



i/p 



{C^K {C{2p/d) - 1 - l/2^p/d) } 



i/p 



C + K- 


X/n=3 (n) 




C + K- 


cm-i-G)^; 



< 1 



□ 

Proposition 1. Let x^y G be positive real values such that y > x. Then 
\/A,B eR, A,B >0 

A^Bx X 
A^By ^ y' 

Proof. Since A^x^y > 0, we can find two positive real values ki > and /c2 > 
such that A = kix and A = k2y. Furthermore, y > x simply implies ki > /c2. 
Using ki and k2 now we can write 

A + Bx kix + Bx ^ /c2X + Bx x 



A^By k2y^By k2y ^ By y 



□ 



References 

[1] W. F. Ames. Numerical Methods for Partial Differential Equations. Aca- 
demic Press, 1977. 

[2] W. E. Arnoldi. The principle of minimized iterations in the solution of the 
matrix eigenvalue problem. Quarterly of Applied Mathematics^ 9(17): 17-29, 
1954. 

[3] E. Bernardis, E. Konukoglu, Y. Ou, D. N. Metaxas, B. Desjardins, and 
K. M. Pohl. Temporal shape analysis via the spectral signature. In In 
Proceedings of MICCAI, 2012. 

[4] F. L. Bookstein, P. D. Sampson, A. P. Streissguth, and P. D. Connor. 
Geometric morphometries of corpus callosum and subcortical structures in 
the fetal-alcohol-affected brain. Teratology, 64:4-32, Jul 2001. 



29 



[5] A. M. Bronstein, M.M. Bronstein, M. Mahmoudi, R. Kimmel, and 
G. Sapiro. A Gromov-Hausdorff framework with diffusion geometry for 
topologically-robust non-rigid shape matching. International Journal of 
Computer Vision, 89(2,3) :266-286, 2010. 

A.M. Bronstein, M.M. Bronstein, and R. Kimmel. Numerical geometry of 
non-rigid shapes. Springer, 2008. 

A.M. Bronstein, A.M. Bruckstein M.M. Bronstein, and R. Kimmel. Analy- 
sis of two-dimensional non-rigid shapes. International Journal of Computer 
Vision, 78(l):67-77, 2008. 

M. M. Bronstein and A. M. Bronstein. Shape recognition with spectral 
distances. IEEE TPAMI, 33(5):1065-71, 2011. 

M.M. Bronstein. Scale-invariant heat kernel signatures for non-rigid shape 
recognition. In Computer Vision and Pattern Recognition, 2010. 

Q.M. Cheng and H.C. Yang. Bounds on eigenvalues of the Dirichlet Lapla- 
cian. Mathematische Annalen, 337:159-175, 2007. 

R. Courant and D. Hilbert. Method of Mathematical Physics, vol I Inter- 
science Publishers, 1966. 

L.C. Evans. Partial Differential Equations. American Mathematical Soci- 
ety, 1998. 

S. Gnutzmann, P. Karageorge, and U. Smilansky. Can one count the shape 
of a drum? Physical Review Letters, 97(9), 2006. 

S. Gnutzmann, U. Smilansky, and N. Sondergaard. Resolving isospectral 
"drums" by counting nodal domains. Journal of Physics A: Mathematical 
and General, 38(41), 2005. 

C. Gordon, D. Webb, and S. Wolpert. Isospectral plane domains and sur- 
faces via riemannian orbifolds. Inventiones Mathematicae, 110(1): 1-22, 
1992. 

P.F.U. Gotardo, K.L. Boyer, J. Saltz, and S.V. Raman. A new deformable 
model for boundary tracking in cardiac MRI and its application to the 
detection of intra- ventricular dyssynchrony. In IEEE Computer Vision and 
Pattern Recognition, 2006. 

N. Iyer, S. Jayanti, K. Lou, Y. Kalyanaraman, and K. Ramani. Three- 
dimensional shape searching: state-of-the-art review and future trends. 
Computer-Aided Design, 37:509-30, 2005. 

V. Jain and H. Zhang. Robust 3d shape correspondence in the spectral 
domain. In International Conference Shape Modeling and Applications, 
2006. 

G. Jurman, R. Visintainer, and C. Furlanello. An introduction to spectral 
distances in networks (extended version). Preprint in ArXiv, Oct 2010. 



30 



[20] M. Kac. Can one hear the shape of a drum? The American Mathematical 
Monthly, 73(7):l-23, 1966. 

[21] R. Lai, Y. Shi, 1. Dinov, T. Chan, and A. Toga. Laplace-beltrami nodal 
counts: A new signature for 3d shape analysis. In International Symposium 
on Biomedical Imaging, 2009. 

[22] R. Lai, Y. Shi, K. Scheibel, S. Fears, R. Woods, A. Toga, and T. Chan. 
Metric-induced optimal embedding for intrinsic 3d shape analysis. In Com- 
puter Vision and Pattern Recognition, 2010. 

[23] B. Levy. Laplace-beltrami eigenfunctions towards an algorithm that "un- 
derstands" geometry. In Shape Modeling and Applications, 2006. 

[24] P. Li and S.T. Yau. On the Schrodinger equation and the eigenvalue prob- 
lem. Communications in Mathematical Physics, 88:309-18, 1983. 

[25] Z. Lian, A. Godil, B. Bustos, M. Daoudi, J. Hermans, S. Kawamura, 
Y. Kurita, G. Lavoue, H. V. Nguyen, R. Ohbuchi, Y. Ohkita, Y. Ohishi, 
F. Porikli, M. Renter, I. Sipiran, D. Smeets, P. Suetens, H. Tabia, and 
D. Vandermeulen. SHREC'll track: Shape retrieval on non-rigid 3D wa- 
tertight meshes. In Proceedings of the Eurographics/ ACM SIGGRAPH 
Symposium on 3D Object Retrieval, 2011. 

[26] T. Mansi, I. Voigt, B. Leonardi, X. Pennec, S. Durrleman, M. Sermesant, 
H. Delingette, A.M. Taylor, Y. Boudjemline, G. Pongiglione, and N. Ay- 
ache. A statistical model for quantification and prediction of cardiac re- 
modeling: application to tetralogy of fallot. IEEE Transactions on Medical 
Imaging, 30(9):1605-16, 2011. 

[27] H.P. McKean and I.M. Singer. Curvature and the Eigenvalues of the Lapla- 
cian. Journal of Differential Geometry, 1:43-69, 1967. 

[28] F. Memoli. A spectral notion of gromov-wasserstein distance and related 
methods. Applied and Computational Harmonic Analysis, 30(3):363-401, 
2010. 

[29] M. Niethammer, M. Renter, F.-E. Wolter, S. Bouix, N. Peinecke, M.-S. 
Koo, and M. E. Shenton. Global medical shape analysis using the Laplace- 
Beltrami spectrum. In Medical Image Computing and Computer Assisted 
Intervention, 2007. 

[30] M. Ovsjanikov, A. M. Bronstein, M. M. Bronstein, and L. J. Guibas. Shape 
google: a computer vision approach to invariant shape retrieval. In Inter- 
national Conference on Computer Vision Workshops, pages 320-327, 2009. 

[31] A. Pleijel. A study of certain green's functions with applications in the 
theory of vibrating membranes. Arkiv For Matematik, 2(6):553-569, 1954. 

[32] M.H. Protter. Can one hear the shape of a drum? Revisited. SIAM Review, 
29(2):185 - 197, June 1987. 

[33] M. Renter, F.-E. Wolter, and N. Peinecke. Laplace-Beltrami spectra as 
'Shape-DNA' of surfaces and solids. Computer- Aided Design, 38:342-66, 
2006. 



31 



[34] M. Reuter, F.-E. Wolter, M. Shenton, and M. Niethammer. Laplace- 
beltrami eigenvalues and topological features of eigenfunctions for statisti- 
cal shape analysis. Computer- Aided Design^ 41:739-55, 2009. 

[35] R. M. Rustamov. Laplace-beltrami eigenfunctions for deformation invariant 
shape representations. In Eurographics Symposium on Geometry Process- 
ing, 2007. 

[36] D.W. Shattuck, M. Mirza, V. Adisetiyo, C. Hojatkashani, G. Salamon, K.L. 
Narr, R.A. Poldrack, R.M. Bilder, and A.W. Toga. Construction of a 3d 
probabilistic atlas of human cortical structures. Neurolmage, 39(3): 1064- 
80, 2008. 

[37] L. Smith. The asymptotics of the heat equation for a boundary value 
problem. Inventiones Mathematicae, 63:467-493, 1981. 

[38] J. Sun, M. Ovsjanikov, and L. Guibas. A concise and provably informative 
multi-scale signature based on heat diffusion. In Eurographics Symposium 
on Geometry Processing, 2009. 

[39] J.B. Tenenbaum, V. de Silva, and J.C. Langford. A global geometric frame- 
work for nonlinear dimensionality reduction. Science, 290(5500) :2319-23, 
2000. 

[40] D.V. Vassilevich. Heat kernel expansion: user's manual. Physics Reports, 
388(5-6) :279-360, 2003. 

[41] H. Weyl. Das asymptotische verteilungsgesetz der eigenwerte linearer par- 
tieller different ialgleichungen. Mathematische Annalen, pages 441-69, 1912. 

[42] E.T. Whittaker and G.N. Watson. A course of modern analysis. Cambridge 
Mathematical Library, 1996. 

[43] H.C. Yang. Estimates of the difference between consecutive eigenvalues. 
preprint, 1991. 

[44] D. Zhang and G. Lu. Review of shape representation and description tech- 
niques. Pattern Recognition, 37(1): 1-19, 2004. 



32 





NN 


FT 


ST 


E 


DCG 


WESD (p^) 
p=3.15, N=100 


0.9933 


0.9020 


0.9305 


0.6900 


0.9706 


WESD (p^^) 
p=2.0, N=100 


0.9933 


0.8923 


0.9238 


0.6824 


0.9691 


N=12, norml 


0.9967 


0.8896 


0.9521 


0.6959 


0.9748 


PSD 
N=12, normA 


0.9917 


0.9153 


0.9569 


0.7047 


0.9783 


N=15, norml 


0.9933 


0.8683 


0.9431 


0.6895 


0.9705 



(a) 




(b) p 3.15 and N 100 



1 

0.99 
0.98 
0.97 
0.96 



-NN 

hDCG 



10 38 66 94 122 
N : number of modes 



1 

0.99 
0.98 
0.97 
0.96 



-»-NN 
-■-DCG 



1.65 2.07 2.49 2.91 3.33 3.75 
p : norm type 




10 38 66 94 122 150 1.65 2.15 2.65 3.15 3.65 



N : number of modes p ; norm type 




38 66 94 122 150 
N : number of modes 



0.695 
0.69 

0.685 
0,68 

0.675 



2.15 2.65 3.15 
p : norm type 



(c) p = 3.15 



(d) N = 100 



Figure 6: Shape-based Object Retrieval Results on SHREC Dataset. a) Re- 
trieval scores obtained by WESD for two different sets of N and p values 
along with the scores obtained by the distance proposed in 33 (values taken 
from 25 .) b) Precision- Recall curves obtained for shape retrieval via WESD 
for the entire dataset. c) Effect of the signature size N on the retrieval scores 
obtained by WESD for a fixed p = 3.15. d) Effect of the norm type p on the 
same scores for a fixed N = 100. 
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(c) Four Putamen (f) p^^^ - surface (g) WESD - surface 

smoothing smoothing 



Figure 7: 2D embedding of subcortical structures: 240 structures (80 caudate 
nucleus, 80 putamen and 80 hippocampus) are extracted from MR scans of 40 
different individuals. (a),(b) and (c) show some example structures from this 
dataset. Note the high intra-class variability and the artefacts due to finite 
resolution and manual segmentations, (d) and (e) plot 2D embeddings of these 
240 structures obtained based on the affinity matrices computed via and 
WESD respectively. These embeddings are computed without any preprocess- 
ing applied to the structures. The embedding obtained with WESD distinctly 
clusters the objects with respect to the anatomical structures. The embedding 
in (d) however, shows some ambiguities in the separation. Graphs in (f) and (g) 
plot the similar embeddings obtained after smoothing the surfaces of the struc- 
tures to remove artefacts. The embedding obtained by Psd^ although better 
than (d), still suffer from similar problems. The embedding based on WESD on 
the other hand, now even between better separates the groups. 
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Figure 8: Analysing 3D + time (4D) cardiac images: Top row shows corre- 
sponding 2D slices of a 4D MRI dataset at time points t = {0, 3, 6, 9, 12}. In the 
middle row, 3D shapes extracted at each of the time points. For five patients, 
we compute the nWESD shape dissimilarity score of the LV blood pool at each 
time point with respect to its shape at t = 0. The graph plots these scores. We 
note that the proposed shape distance is able to capture the dynamic process of 
the LV shape changes and furthermore, the symmetry between the two phases 
of an heart beat: diastole and systole. 
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